setwd("D:/PaperRev2/Plot")
#This is a script for R to calculate column density from Planck Observational Data, following A&A_536_A23P.#
#---------The Constants-----------------#
h<-6.63e-27        #(erg.s)
kB<-1.38e-16       #(erg.k^{-1})
c<-2.998e10        #(cm/s)
nu<-8.57e11        #(Hz)
kappa_std<-0.1     #At 1 THz          
nu_std<-1e12       #(Hz)
mJy<-1e-23*1e-3    #(erg/(s cm^2 Hz)) 
mu<-2.33           #Mean weight of molecules
mH<-1.67e-24       #(g), Weight of proton
arcmin<-1/60/180*pi#arcmin->rad
#---------------------------------------#
#---------The Functions-----------------#
PlanckLaw<-function(freq,t){
2*h*freq^3/(c^2)/(exp(h*freq/kB/t)-1)
}
Opacity<-function(beta,nu,kappa_std,nu_std){
kappa_std*((nu/nu_std)^beta)
}        
#---------------------------------------#
#---------Input the Data----------------#  
rt<-read.table("ECC_All.txt")
kappa<-Opacity(2,nu,kappa_std,nu_std)
name<-rt$NAME
flux<-rt$APFLUX857
flux<-flux*mJy
T<-rt$TEMPERATURE_CORE
B<-PlanckLaw(nu,T)

smax<-rt$MAJ_AXIS_FWHM_CORE 
#smax<-4.42/2.355

smax<-smax*arcmin/2.355

smin<-rt$MAJ_AXIS_FWHM_CORE 
#smin<-4.42/2.355

smin<-smin*arcmin/2.355

Omega<-pi*smax*smin
beta<-rt$BETA_CORE
NH2<-flux/(Omega*kappa*B*mu*mH)


#This is a script for R to calculate column density from PMO CO isotopes data.#

list<-read.table("TPC_list.TXT", header=TRUE)
dust<-read.table("ECC_All.TXT")
i<-1
while(i<(length(list[[1]])+1))
{
		j<-1
		name<-as.character(list[i,1])
		date<-as.character(list[i,2])
		lname<-as.character(list[i,3])
		dir<-"D:/TaurusBUR/Stat"
		dir<-paste(dir,name,sep="/")
		dir<-paste(dir,date,sep="")
		dir<-paste(dir,"NH2-out_1.txt",sep="/")
    gas<-read.table(dir,header=FALSE)
    ra<-gas[[1]]
    dec<-gas[[2]]
    pNH2<-gas[[3]]/1e21
    while(lname!=as.character(dust[j,1])&j<(length(dust[[1]])+1))												
    {	
    	j<-j+1									
    }	
    dust_name<-as.character(dust[j,1])											
    dust_max<-dust$MAJ_AXIS_FWHM_CORE[j]*60/2.355
    dust_min<-dust$MIN_AXIS_FWHM_CORE[j]*60/2.355
    dust_r2<-dust_max^2+dust_min^2
    j<-1
    l<-1
    gas_NH2<-0
    while(j<(length(ra)+1))
    {
    	if((ra[j]^2+dec[j]^2)<dust_r2)
    		{
    			gas_NH2<-gas_NH2+pNH2[j]
    			l=l+1
    		}
    	j<-j+1
	  }
	  gas_NH2<-gas_NH2/l
	  list[i,4]<-lname
	  list[i,5]<-gas_NH2
	  list[i,6]<-NH2[j]/1e21
	  list[i,7]<-l
	  i<-i+1
}
SourceName<-list[[4]]
NH2_Gas<-list[[5]]
NH2_Dust<-list[[6]]
TPC_out<-data.frame(SourceName,NH2_Gas, NH2_Dust)
write.table(TPC_out,file="TPC_out.dat")


